The Surface Energy of a Bounded Electron Gas: Analysis of the Accuracy of the 
Local-Density Approximation via Ab Initio Self-Consistent-Field Calculations 
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We report an ab initio evaluation of the surface energy 
of a simple metal, performed via a coupling-constant integra- 
tion over the dynamical density-response function. The rapid 
rate of change of the electron density at the surface is treated 
exactly. Long-range correlations are treated self-consistently 
in the random-phase approximation; short range correlations 
are included in time-dependent local density-functional the- 
ory. Our results provide a numerical measure of the error 
introduced by the usual local-density approximation; this er- 
ror is found to be small. 



Since the pioneering work of Lang and KohnEJ, the cal- 
culation of the surface energy of a metal has been the 
subject of long-standing interest. These authors were the 
first to include the crucial effects of exchange and correla- 
tion self-consistently within the local-density approxima- 
tion (LDA) of density-functional theory (DFT)H. Lang 
and Kohn also discussed the effect of the crystal lattice, 
whose full inclusion within the LDA poses no difhcul- 
ties of principle these days. By contrast, the question of 
the impact of non-local Coulomb correlations, and their 
interplay with the strong charge inhomogeneiijj^ at the 
surface, has remained unsettled over the yearaj^lj. 

Recent calculations have rekindled the controversy on 
the question of the quality of the LDA surface energy. 
A many-body scheme which starts from a physically- 
motivated model of the correlated ground state, and 
treats the interactions "wia a Fermi hypernetted-chain ap- 
proximation (FHNCA)t3, has yielded surface energies 
which are signifi|Cantly higher than the LDA results of 
Lang and KohrO. By contrast, the surface energies ob- 
tained in density-functional calculationsllj baaed on the 
use of the Langreth-Mehl non-local functional^ are much 
closer to the LDA result. Finally, althougt. very recent 
quantum Monte Carlo (QMC)calculationf£j ptffree with 
the latter non-local density functional resultstHi for high 
densities (r^ < 2.07), they agree with the FHNCA for 
lower densities (r^ > 3.25). 

The purpose of this paper is to establish, in a controlled 
way, the impact of non-locality on the surface energy of 
an electron gas. To this end we carry out non-local and 
local calculations within exactly the same conditions, i.e., 
we consider the LDA as a special case of the general non- 
local formalism based on a coupling-constant integration 



over the dynamical density-response function. First, the 
effects of non-local correlations are investigated fully self- 
consistently within a well-defined many-bady framework, 
the random-phase approximation (RPA)lli. Our side-by- 
side calculation, in which the same diagram is used to 
generate the local and non-local surface energies, shows 
that the LDA is quite accurate over the entire density do- 
main appropriate to metals (rs = 2 — 6). Of course, the 
absolute values of our RPA energies cannot be expected 
to be more reliable than, say, the QMC surface energies. 
However, the significance of our results lies in the elucida- 
tion of the difference between non-local and local surface 
energies. We also explore the impact of short-range cor- 
relations by invoking a time-dependpHt extension of local 
density-functional theory (TDLDA)^. The TDLDA ver- 
tex introduces an element of arbitrariness, since it con- 
tains no dynamical effects. Thus, again, we are less inter- 
ested in the absolute value of the TDLDA surface energy 
than we are in its difference with its local counterpart; 
again, we evaluate this difference in an unambiguous way. 
Our TDLDA results support the conclusion drawn from 
our RPA results that the error introduced by the LDA is 
not large. 

The ground-state energy of an interacting electron sys- 
tem can be witten as a functional of the electron number 
density ri(r)a. 



E[n] = Ek[n] + E^^H + E^c[n], 



(1) 



where -Efe[n] is the kinetic energy of a non- interacting 
system with the same density and £'es['T-] is the Hartree 
electrostatic energy. The exchange-correlation (XC) en- 
ergy, Exc[n], can be obtained from a coupling constant- 
integration over the interaction energy. We follow Lan- 
greth and PerdewD and choose the coupling constant such 
that the density is maintained at its fully-interacting 
value while the electron-electron interaction strength is 
varied from A = to A = 1. The fluctuation-dissipation 
theorem then leads us to the result that 
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where x;^(r, r';w) is the density-response function. 

If the Coulomb correlations are ignored altogether, Eq. 
(2) reduces to the expression 
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is the density-response function for non-interacting elec- 
trons. For inhomogeneous systems Eq. (3) coincides with 
the Fock exchange energy only if the wave functions V'i (r) 
are the solutions of the non-local Hartree-Fock equation. 

Actually, Eq. (4) can be interpreted in a more gen- 
eral context as giving the density-response to an ap- 
propriate mean field set up by the dynamical polar- 
ization of the Fermi sea. In the particular case of 
the RPA, the single-particle wave functions ipi{r) en- 
tering Eq. (4) are strictly the self-consistent eigen- 
functions of the one-electron Hartree Hapailtonian. In 
time-dependent density- functional theoryll3, the "non- 
interacting" electrons in question are described by the 
solutions of the time-dependent counterpart of the Kohn- 
Sham equationcl; in usual practice, these amplitudes 
are approximated by standard LDA wave functionslHI. 
Both approaches to the evaluation of the polarizability 
x''(r, r';w) will be considered in our numerical study of 
the surface energy of a simple metal. 

In both RPA and TDLDA the response function satis- 
fies the integral equationllj 

XA(r,r';w) =x°(r,r';w) 

-h / dri / dr2X°{r,ri;uj)XV{ri,r2)x\{r2,-r';uj). (5) 

In the RPA, the effective electron-electron interaction 
F(ri, r2) is just the bare Coulomb potential. In TDLDA, 
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Vxc being the local XC potential. The combination of 
Eqs. (2) and (5) defines either the TDLDA or the RPA 
exchange-correlation energy, depending on whether or 
not the Coulomb interaction is "dressed" according to 
Eq. (6). This dressing corresponds to the inclusion of 
short-range correlations, which are ignored in RPA. 

We consider a jellium slab of thickness a and density 
77+ — qp/Sn'^, where qp — {9n/4)^/^/{rsao) is the Fermi 
wave vector, r^ is the Wigner-Seitz radius, and oq is the 
Bohr radius. The slab is translationally invariant in the 
plane of the surface, which is assumed to be normal to 
the z axis. Thus the single-particle wave functions are of 
the form 






(7) 



where qy is a wave vector parallel to the surface and 
A is the normalization area. The wave functions 4>i{z) 
describe motion normal to the surface, and are obtained 
self-consistently with the effective one-electron potential, 
Veff{z). In the RPA, Vejf consists of just the Hartree 
potential Ves{z); in TDLDA, 



Veff{z) = Ves{z)+VUz). 



(8) 



For reference, we recall that, in the simplest non-self- 
consistent microscopic model of the surface, the infinite- 
barrier model (IBM), Veff{z) is replaced by an infinite 
square barrier, and the functions (j)i{z) are simply sines. 
A solution of Eq. (5) for a self-consistent description 
of the surface potential has been given some time agcO. 
We assume that n{z) vanishes at a distance zq from ei- 
ther jellium edge, and expand the wave functions 4>i{z) 
in a Fourier sine series; zq is chosen sufficiently large for 
the physical results to be insensitive to the precise value 
employed. We introduce a double-cosine Fourier repre- 
sentation for the density response function, and also for 
the Coulomb potential and the Dirac delta function en- 
tering Eq. (2)E3. The use of this representation allows us 
to perform analytically the integrals in Eq. (2) involv- 
ing the coordinate normal to the surface; the integrals 
over parallel-momentum transfers, over energy transfers, 
and over the coupling constant, are performed numeri- 
cally. The total energy given by Eq. (1) is evaluated 
in a similar way. Subtracting from the total energy the 
corresponding result for a homogeneous electron gas of 
density n+, E^ \n+], we obtain the surface energy 
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The LDA is obtained from the above non-local formal- 
ism by replacing the response function entering Eq. (2) 
by its counterpart for a homogeneous electron gas with 
the local value ofr-the density. This replacement leads us 
to the result thatc3 
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Equation (10), with E^^ [n{z)] evaluated on the basis of 
(homogeneous electron gas-) RPA and TDLDA density- 
response functions, calculated for the local value of the 
density n{z) obtained self-consistently with Hartree and 
LDA effective potentials, yields our RPA- and TDLDA- 
based LDA surface energies. We will consider these local 
results together with the LDA results of Lang and Kohntl 
—which we also obtain from Eq. (10) through the use of 
the Wigner interpolation formula for E^^ [n{z)'f3. 

Great care was exercised to ensure that our slab calcu- 
lations are a faithful representation of the surface energy 
of a semi-infinite medium. This issue is important, in 
view of the subtle cancellations which exist between the 
various contributions to the surface energy; furthermore, 
these contributions are oscillatory functions of the slab 



width a. (The amphtude of the oscillations decays ap- 
proximately linearly with a; their period equals Ai?/2, 
Xp — 2Tr/qp being the Fermi wavelength.) 

For each value of rg we have actually considered three 
different values of a. One such value, a„, is the thresh- 
old width for which the n-th subband for z-motion is 
first occupied; for this width the surface energy is a lo- 
cal minimum. The other two values of the slab width, 
a~ = an — Xf/4:, and a+ = a„ + Af/4, correspond to 
the twOplocal maxima about the minimum. Utilizing the 
relatione^ 



(11) 



we are able to extrapolate our calculated surface en- 
ergy to the infinite- width limit. This procedure was first 
tested, with values of n up to 200, for the IBM, for which 
analytical insight is possible by virtue o£_thc simple na- 
ture of the one-electron wave functionso. The results 
presented below correspond to slabs with n = 12, for 
which a « 5 — 6Af , depending on Ts . Based on this pro- 
cedure, we estimate that the numerical error introduced 
by our slab simulations corresponds to one unit in the 
last digit of all the entries in Table I. (We remark that 
our results were found to be insensitive to the precise 
value of the number of sines, Smax, kept in the expansion 
of the wave functions 4'i{z), for Smax > 280.) 

The key results of our work can be readily grasped from 
Fig. 1, in which we show the surface energy as a func- 
tion of rs- Consider first the RPA. The reasons for the 
significance of our RPA calculations are: (i) the effects of 
long-range correlations are included fully self-consistently 
with the electron density profile (which, we recall, is eval- 
uated in the Hartree approximation); (ii) the non-local 
and local calculations are carried out within one and the 
same density-response framework; (iii) this framework is 
devoid of any ambiguities in the treatment of the many- 
body problem. It is clear from Fig. 1 that the local RPA 
surface energy differs little from its non-local counterpart 
over the entire metallic range of densities. 

Next, we consider the effects of the short-range correla- 
tions built into the XC potential, V^c- In the full TDLDA 
treatment, this effect is included in both the one-electron 
potential of Eq. (8) and the electron-electron interaction 
of Eq. (6). (We evaluate-yrc with use of the Perdew and 
Zunger parametrizationE3.) Overall, the main impact of 
the inclusion of XC is via Eq. (8), which, through self- 
consistency, yields a more abrupt electron density profile 
at the surface, relative to the Hartree profile. This leads 
to the large lowering of the surface energy, relative to the 
RPA, which we observe in Fig. 1. For completeness, in 
Fig. 1 we show both the full TDLDA (solid line), and 
the result obtained upon including XC in Eq. (8) but not 
in Eq. (6) (dashed line); clearly, the impact of the XC 
vertex is smaller than the effect of the inclusion of XC in 
the electron density profile. 

As was the case above with the RPA calculations, the 
difference between non-local (solid line) and local (dotted 



line) TDLDA surface energies is relatively small. More 
specifically, the error introduced by the LDA is of the 
order of 50% smaller than the error one would impute 
to the local apjpcpximation on the basis of the non-local 
FHNCA resultalj (open circles in Fig. 1), particularly in 
the crucial high-density region (r^ « 2). 

It is apparent from Fig. 1 that our non-local TDLDA 
surface energies agree Wjell, for all densities, with those 
obtaijifid by Zhang et al.ll3 using the non-local Langreth- 
MehO XC functional. By contrast|,— while the surface 
energy obtained via QMC techniqueslij is close to either 
set of results for rg — 2.07, for lower densities the QMC 
surface energies are appreciably larger; in fact, they are 
close to the FHNCA values. 

It is interesting to note that our TDLDA-based LDA 
surface energies agree [Closely with the LDA calculations 
of Acioli and CeperleyHI. (Our results can also be repro- 
duced from Eq. (10) with usG_of the Perdew and Zunger 
parametrization of £';^[n(z)]Ej.) Thus, while from the 
QMC results of Ref. |l^ one would conclude that the error 
introduced by the LDA is, for r^ > 3.25, rather signifi- 
cant, this is not what transpires from our results —most 
particularly, from our RPA calculations. 

We note that the difference between our TDLDA- 
based LDA results and the local Lang-Kohn surface ener- 
gies (dashed-dotted-dotted-dotted lipe) is simply a con- 
sequence of the use foi these authora2jof the Wigner for- 
mula for correlationEll. The much larger difference be- 
tween Lang-Kohn surface energies and our RPA-based 
LDA results reflects, again, the impact of XC on the elec- 
tron density profile at the surface. 

If the surface energy is obtained from Eq. (3) (not Eq. 
(2)), in conjunction with an exchange-only correction to 
the Hartree barrier in Eq. (8), we obtain non-local sur- 
face energies which are closcj-to the Hartree-Fock results 
reported by Krotscheck et al.t3 By contrast, the non-local 
correlation contribution to the surface energy reported in 
Ref. no is significantly higher than our ab initio correla- 
tion surface energy; as a result, the FHNCA surface ener- 
gies are much higher than our non-local TDLDA values, 
as shown in Fig. 1. The FHNCA and RPA surface ener- 
gies turn out to be quite close, at low densities, because 
of compensation between the effect of XC in the effec- 
tive one-electron potential, which is absent in the RPA, 
and the very large FHNCA non-local correlation surface 
energy. 

In conclusion, we have presented ab initio calcula- 
tions of the surface energy of a bounded electron gas. 
The unambiguous nature of the comparison of local- 
vs. non-local surface energies made possible by our self- 
consistent RPA calculations, leads us to the conclusion 
that the LDA does, within the RPA, work. We have 
also evaluated the difference between non-local and lo- 
cal TDLDA surface energies, and the results so obtained 
support the conclusion that the error introduced by the 
LDA is, within the TDLDA, not large, over the whole 
range of electron densities appropriate to metals. Fur- 
ther progress in the quantitative ab initio evaluation of 



this difference requires improvements in the treatment of 
dynamical many-body correlations beyond the scope of 
Eq. (6). 
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TABLE L Exchange-correlation [oxc) and total (a) 
non-local surface energies obtained in RPA and TDLDA, and 
their local counterparts. Units are erg/cm? . 
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FIG. 1. Non-local RPA and TDLDA surface energies 
(solid lines), as functions of r^,. The dashed line is the 
TDLDA result obtained upon excluding the XC vertex from 
Eq. (6). Dashed-dotted and dotted lines represent the lo- 
cal RPA and TDLDA surface energies, respectively. The 
dashed-dotted-dotted-dotted line represents the Lang-Kohn 
surface energy of Ref. hi. Open circles, stars and open squares 
are results taken from Refs. hol, hdand [l3L respectively. 
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